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Abstract. GdN bulk is studied with the local density approximation, on the Hartree- 
Fock level, and on the level of the hybrid functional B3LYP. A local basis set formalism 
is used, as implemented in the present CRYSTAL06 release. It is demonstrated that the 
code is technically capable of treating this system with its 4/ electrons explicitly, i.e. 
out of the core. The band structure at the level of the local density approximation is in 
good agreement with earlier calculations and is found to be half-metallic. The Hartree- 
Fock band structure is insulating with a large gap. Interestingly, three solutions were 
found at the B3LYP level. The lowest of them is insulating for majority spin, and the 
Fermi surface for minority spin consists only of points, resulting in a very low density 
of states around the Fermi level. 
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1. Introduction 

GdN crystallizes in a rock salt structure and is by now established as a ferromagnet with 
a Curie temperature of ~ 69 K [H [2j [31 H] . The question whether it is an insulator or 
not is still under discussion. In the transmittance spectrum, a gap of 1 eV was observed 
[5]. From optical reflection and derivation of the plasma resonance and from the Hall 
effect |6] it was argued that GdN was a semimetal. Thin GdN layers were found to be 
insulating in resistivity measurements (?]. In X-ray photoelectron spectroscopy (XPS) 
experiments, the occupied 4/ bands were found at a binding energy of 7.8 eV for GdN 
[8]. For the related systems GdX (X=P, As, Sb, Bi), the occupied Gd 4/ states were 
found at ~ 9 eV, and the unoccupied states at ~ 5 eV in XPS and X-ray bremsstrahlung 
isochromat spectroscopy (X-BIS) measurements [9]. 

The first theoretical study of GdN was an augmented plane wave calculation using 
the Slater exchange potential [ID], which resulted in an insulator with a very small gap 
[II] . The / electrons were treated as core states in this calculation. A further calculation 
with the / electrons in core [12] using the local density approximation (LDA) gave a 
similar band structure, but without gap. The spin-polarized solution was found to be 
metallic for majority, and insulating for minority spin. Recently, several calculations 
with explicitly treating the / electrons were performed. When the pure LDA was used, 
a metal was obtained [131 [T4] , with the occupied / states at ~ -4 to -3 eV below the 
Fermi level. To take into account the strongly correlated nature of the / electrons, a 
Hubbard Uf term was applied to these states [15j [El [H] . This still resulted in a metallic 
solution, and only when also a Ud term was applied to the Gd d states, an insulator 
was obtained [13] . Similarly, in [15], LDA+U calculations gave a metallic ground state, 
and by applying further rigid shifts to the 5d and 4/ states, an insulating state was 
found. When the lattice was expanded, also a transition to an insulator was predicted, 
at the LDA+U level [16] . A slightly different approach to take into account the strongly 
correlated nature of this system is to apply a self-interaction correction. This led again 
to a half- metallic ground state for GdN, with a gap in the minority channel [17] . 

Therefore, there are many reasons to study GdN. Besides these aforementioned 
reasons, from the theory point of view, a further point is that it would be interesting 
to test the performance of hybrid functionals, as a further method to treat strongly 
correlated systems. It has been shown that they give surprisingly good values for the 
band gaps for such systems, where the local density approximation and other standard 
functionals often fail [HI [191 120] . Recently, the B3LYP hybrid functional was applied to 
UO2, with explicitly treated / electrons, and also a good value for the gap was obtained 
[21] . Also, plutonium oxides were found to be well described by hybrid functionals [22J. 
Cerium oxides were studied using screened hybrid density functionals with a local basis 
set [23], and with a plane wave basis set and hybrid functionals [24J. 

The CRYSTAL code used in the present work is based on a local Gaussian basis set. 
The first version released in 1988 was a pure Hartree-Fock code. Hybrid functionals, 
which use an admixture of exact (Fock) exchange, were available from the 1998 release 
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onwards. It turned out that, in the CRYSTAL implementation, the CPU times for 
calculations with hybrid functionals are comparable to the CPU times for calculations 
with standard functionals such as the LDA or gradient corrected functionals, and also 
the memory requirements are similar in both cases. 

The present CRYSTAL release (25] was announced of being able to use /-functions 
as polarization functions. In this article, it will be shown that it is also technically 
possible to perform calculations on systems with / electrons, with GdN as an example. 
There are various reasons why GdN was chosen as a test case: firstly, the /-occupancy 
if / 7 , and the half-filled /-shell is probably one of the easiest examples to start with 
calculations on /-electron systems. Also, the crystal structure is fairly simple (NaCl 
type). In addition, as mentioned earlier, the system has already been theoretically 
described with other codes. This gives the opportunity to compare the present 
calculations, where possible, with the results obtained earlier. 

2. Method and calculational details 

The calculations were performed with the CRYSTAL06 code [25] which uses a local 
basis set. GdN was considered in the face centered cubic lattice. A ferromagnetic order 
was assumed, as this is the experimentally observed state. 

The local density approximation, the Hartree-Fock approach (HF), and the B3LYP 
hybrid functional were employed. The local basis sets were chosen as follows: for 
Gd, a small core pseudopotential [26] was used, which includes scalar-relativistic 
effects (mass- velocity, Darwin and averaged spin-orbit operator [27]). Together with 
the pseudopotential, the corresponding basis set was used [28], with the following 
modification: the inner [8s7p4d4f] were kept as in the original basis set, and in addition 
one diffuse s, p and d function with exponent 0.12 were added, so that a [9s8p5dAf] basis 
set was obtained. The [3s2pld] nitrogen basis set from reference [29] was used (6-21G 
as in [30] and a d-function with exponent 0.8). The present release of CRYSTAL06 
can not compute the atomic solution of atoms with occupied /-orbitals. This solution 
would normally be used as an initial guess for the subsequent calculations on the 
periodic system. To overcome this problem, the atomic occupancy of Gd was chosen 
as 4s 2 4p 6 4<i 10 5s 2 5p 6 5(i 8 6s 2 which is obviously unusual, but sufficient to overcome the 
problem of having /-orbitals occupied (ls 2 2s 2 2p 6 3s 2 3p 6 3c/ 10 electrons are described by 
the pseudopotential). For nitrogen, the ls 2 2s 2 2p 3 state was occupied. No convergence 
is achieved for the Gd atom which is not too surprising due to the unusual occupancy, 
but this does not cause severe problems either. To converge GdN bulk, one convergence 
strategy is then e.g. to use the Anderson mixing scheme [31] with 80% mixing, and in 
addition to keep the spin locked to a value of 7/2 in the first few (e.g. 10) cycles. In 
some cases, slight variations of this approach were necessary, but in the end convergence 
could be achieved in all cases. The cases of LDA and HF are more straightforward, 
whereas B3LYP has the additional difficulty that three different solutions were found. 
To compute the corresponding potential curves and determine the minimum, one certain 



Electronic structure of GdN, and the influence of exact exchange 4 

Table 1. The computed equilibrium lattice constant and bulk modulus of GdN, at 
various levels of theory. 





a 


B 




(A) 


(GPa) 


HF 


5.10 


167 


LDA 


4.91 


174 


LDA (ref. [12]), 4/ in core 


4.999 


188 


LDA (ref. ,13]), 4/ explicit 


4.98 




B3LYP (solution with lowest energy) 


5.10 


137 


exp. 


4.99 [34] 


192 ± 35 [35] 



solution was used as the initial guess for the other geometries. This way, the potential 
curves on the B3LYP level could be computed for all three solutions found. 

A sampling with 16 x 16 x 16 ^-points in the reciprocal lattice was used. In the 
case of B3LYP, calculations without smearing and with a smearing of 0.01 Eh {Eh are 
hartree units, 1 .£^=27.2114 eV) were performed, to explore possible changes due to 
the different decay properties of the density matrix in metallic systems [32] . However, 
no severe differences between the calculations with and without smearing temperature 
were observed. 

The remaining parameters such as grids for the numerical density functional 
integration were chosen as the default. The only parameter which needed to be varied 
were the various thresholds for integral selection [33J . 

3. Results 

3.1. LDA calculations 

The computed LDA equilibrium lattice constant of 4.91 A and the bulk modulus of 174 
GPa (table CD) agree reasonably well with the literature. The LDA band structure is 
displayed in figure [H It is half-metallic, because the bands for majority spin cross the 
Fermi energy, whereas the minority bands have a gap. The /-bands are positioned at 
~ -3.1 eV below and ~ 2.6 eV above the Fermi energy. The remaining bands are as 
follows: the lowest band displayed is N s, followed by three N p bands which hybridize 
with the Gd d bands which are next higher in energy. The band structure is in good 
agreement with the recent linearized muffin tin orbital approach [131 EH] • 

The Mulliken population analysis gives a charge transfer of -0.7 to the nitrogen 
atoms. The Gd / population is 7.1, the 5d population is 1.3, and the 6s and 6p 
populations are 0.4 each. The spin population is 7.1 on Gd (6.9 due to the / electrons 
and 0.2 due to the d electrons) and -0.1 on N (due to the p states), which agrees well 
with earlier LDA calculation where the /-states were treated in core |12j . 

The corresponding density of states is displayed in figure [2j In the selected energy 
range, the mainly occupied states are N (s and p), Gd d and Gd /. The density of 
states, projected on these states is also shown. 
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Figure 1. LDA band structure of GdN at the computed equilibrium lattice constant 
of 4.91 A. The Fermi energy is positioned at 0. 

3.2. HF calculations 

The Hartree-Fock equilibrium lattice constant is 5.10 A and thus slightly larger than 
the experimental value, which was also observed in an earlier Hartree-Fock calculation 
where the 4/ electrons were treated as core electrons |36j. Similarly, the bulk modulus is 
slightly too small due to this overestimation of the lattice constant. The band structure 
at the Hartree-Fock level is displayed in figure [3j GdN is insulating on the HF level, with 
a gap of ~ 5 eV. The occupied / bands are pushed downwards to ~ — 15 eV below the 
top of the valence bands, and similarly the unoccupied / bands are pushed upwards to 
~ 22 eV. This is to be expected, since any screening, which is present on the LDA level, 
is missing at the HF level, and thus the localized / electrons feel the full, unscreened 
Coulomb repulsion. Another striking feature is that the N p bands and the Gd d bands 
are now well separated. This is also visible in the corresponding density of states in 
figure HI 

The total Mulliken population of N is -1.1 |e|. The spin population is 7.2 on Gd 
(7.0 due to the Gd / electrons, and the remaining 0.2 mainly due to the Gd d states), 
and -0.2 for N (practically completely due to the N p states). 
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Figure 2. LDA density of states at the computed equilibrium lattice constant of 
4.91 A. Besides the total density of states, the projected density of states is shown for 
projections on N, Gd d, Gd f states. The Fermi energy is positioned at 0. The upper 
part of each panel is for majority spin, the lower part for minority spin. 

4. B3LYP calculations 

The B3LYP calculations provided the somewhat surprising result that two states were 
found, which are close in energy. The state with the lowest energy will be discussed first. 
The computed lattice constant of this state is 5.10 A. Its band structure is displayed in 
figure El and the corresponding density of states in figure El The / bands are positioned 
at ~ -5 eV below and ~ 7 eV above the Fermi energy, which agrees better with the 
experiment, compared to LDA and HF. In this case, the majority bands have a gap, in 
contrast to the LDA. The minority valence and conduction bands only touch in various 
places, so that the density of states at the top of the valence band practically vanishes. 

The N population is -0.8 |e|. The Gd spin population is 6.8 (7.0 due to the / 
electrons, 0.1 due to the d electrons and -0.3 due to the Gd s basis functions), and the 
N spin population is 0.2. 

The second state is 0.004 Eh (0.1 eV) higher in energy, with a minimum at 5.11 
A, and a bulk modulus of 133 GPa. Its band structure is displayed in figure 0, and the 
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Figure 3. HF band structure of GdN at the computed equilibrium lattice constant of 
5.10 A. The top of the valence band is positioned at 0. 

corresponding density of states in figure El The / bands are similarly to those of the 
first state (the one lowest in energy) positioned at ~ -5 eV below and ~ 7 eV above 
the Fermi level. This state is half-metallic, in the same way as the LDA solution: the 
majority bands do not have a gap, whereas the minority bands have a gap. 

The N population is also -0.8 |e|. The Gd spin population is 7.4 (7.0 due to the 
/ electrons, 0.1 due to the d electron and 0.3 due to the Gd s basis functions). The N 
spin is -0.4 and thus antiparallel to the total Gd spin, as in the LDA. 

The fact that there are two solutions close is energy is apparently due to the large 
overlap of the / bands with the mainly nitrogen bands, especially at the L point where 
the nitrogen band minimum coincides with the Gd / bands. This overlap also explains 
why the spin populations are a bit different from the LDA and HF values. Having the 
nitrogen spin parallel or antiparallel does apparently not make a huge energy difference 
at the B3LYP level, and these states are close in energy. 

Even a third solution is found, with its band structure displayed in figure [9j This 
solution is however 0.089 Eh (~ 2.4 eV) higher in energy than the lowest solution 
(the equilibrium lattice constant would be slightly shorter, 5.05 A, which does however 
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Figure 4. HF density of states at the computed equilibrium lattice constant of 5.10 
A. Besides the total density of states, the projected density of states is shown for 
projections on N, Gd d, Gd f states. The top of the valence band is positioned at 0. 



not have a huge impact, and the properties of this state at 5.10 A are similar). The 
computed bulk modulus of this solution is 159 GPa. Still, it is interesting to study this 
solution: the / bands are now separated from the nitrogen bands, and this solution is 
insulating, with a gap of ~ 0.7 eV. The charge transfer to nitrogen is ~ -0.9, and the 
spin is now more similar to the LDA and HF solutions: 7.1 for Gd in total (7.0 due 
to the / electrons, 0.1 due to the d electrons), and -0.1 for N due to the p electrons. 
This confirms that the somewhat unusual population for the two lowest states on the 
B3LYP level are due to the hybridization at the L point. This third solution has also 
the property of "interpolating" between LDA and HF (lattice constant, shape of the 
band structure), what is often observed for the B3LYP hybrid functional: due to the 
admixture of Fock exchange, the results for B3LYP calculations are usually expected to 
be in the range in between LDA and HF. 



Electronic structure of GdN, and the influence of exact exchange 9 



GdN: B3LYP, lowest energy solution 

ALPHA BETA 




G X W L G KXG X W L G KX 



Figure 5. B3LYP band structure of GdN at the computed equilibrium lattice constant 
of 5.10 A, for the energetically most favorable solution. The Fermi energy is positioned 
at 0. 

5. Conclusion 

It was demonstrated that Gaussian type orbitals are technically capable for performing 
calculations on GdN bulk with / electrons explicitly treated. The LDA results are in 
very good agreement with previous calculations and give a half-metallic solution, with 
the majority bands being conducting and the minority bands insulating. The Hartree- 
Fock solution is insulating with a large gap of ~ 5 eV, which is a typical overestimation 
of the gap at the Hartree-Fock level due to the lack of screening. On the B3LYP level, 
three solutions were found. The lowest one has a gap for majority spin; and for minority 
spin, the valence and conduction bands only touch at certain points of the Brillouin zone. 
The corresponding density of states is thus very small around the Fermi energy. The 
second solution is very close in energy (~ 0.1 eV higher), and the majority bands cross 
the Fermi energy. A third solution was found to be insulating, but this third solution 
is about 2.4 eV higher in energy than the lowest solution. The fact that there are two 
nearly degenerate solutions makes the comparison with the experiment very difficult, 
and it is difficult to judge the performance of the hybrid functional B3LYP in this case. 
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GdN: B3LYP density of states 
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Figure 6. B3LYP density of states at the computed equilibrium lattice constant of 
5.10 A, for the energetically most favorable solution. Besides the total density of states, 
the projected density of states is shown for projections on N, Gd d, Gd f states. The 
Fermi energy is positioned at 0. 

It is however also very interesting that two different solutions with very different Fermi 
surfaces are obtained in the calculations. As the experimental situation is not fully 
clear, experiments such as photoemission on this system might be very interesting to 
further elucidate the electronic structure of GdN (indeed, spin and angle resolved inverse 
photoemission spectroscopy was also suggested in another recent theoretical work [37]). 
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GdN: B3LYP, second solution 
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Figure 7. B3LYP band structure of GdN at the computed equilibrium lattice constant 
of 5.10 A. This solution is slightly higher (0.1 eV) in energy than the solution in figure 
[5] The Fermi energy is positioned at 0. 
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GdN: B3LYP density of states 
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Figure 8. B3LYP density of states at the computed equilibrium lattice constant of 
5.10 A. This solution is slightly (0.1 eV) higher in energy than the lowest one in figure 
[5] Besides the total density of states, the projected density of states is shown for 
projections on N, Gd d, Gd f states. The Fermi energy is positioned at 0. 
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GdN: B3LYP, third solution 
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Figure 9. B3LYP band structure of GdN at the lattice constant of 5.10 A, for the 
insulating solution. This solution is the highest found at the B3LYP level, about 2.4 
eV higher than the energy of the lowest solution in figure El The middle of the gap is 
positioned at 0. 
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